Self-force with (3+1) codes: a primer for numerical relativists 
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Prescriptions for numerical self-force calculations have traditionally been designed for frequency- 
domain or (1+f) time-domain codes which employ a mode decomposition to facilitate in carrying 
out a delicate regularization scheme. This has prevented self-force analyses from benefiting from the 
powerful suite of tools developed and used by numerical relativists for simulations of the evolution of 
comparable-mass black hole binaries. In this work, we revisit a previously-introduced (3+1) method 
for self-force calculations, and demonstrate its viability by applying it to the test case of a scalar 
charge moving in a circular orbit around a Schwarzschild black hole. Two (3+1) codes originally 
developed for numerical relativity applications were independently employed, and in each we were 
able to compute the two independent components of the self-force and the energy flux correctly 
to within < 1%. We also demonstrate consistency between i-component of the self-force and the 
scalar energy flux. Our results constitute the first successful calculation of a self- force in a (3+1) 
framework, and thus open opportunities for the numerical relativity community in self-force analyses 
and the perturbative modeling of extreme-mass-ratio inspirals. 
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I. INTRODUCTION 



A pressing challenge in gravitational wave source mod- 
eling is the inspiral of a solar mass compact object into 
a supermassive black hole, better known as an extreme- 
mass-ratio inspiral or EMRI. These inspirals result from 
scattering processes in the star-rich cores of galaxies, and 
tend to be highly-eccentric in the strong field region of the 
supermassive black hole [l[. The intricate gravitational 
waves they emit are believed to be the most complicated 
among LISA sources, and their detection and analysis 
promise significant science returns for relativistic astro- 
physics and general relativity 0. For this to come to 
fruition, precise models of their gravitational waves will 
be necessary. 

Immediately confronting this objective are the dramat- 
ically different scales that characterize EMRIs. First, 
there is the short length scale of the distortion on the 
background spacetime made by the compact object, 
which will need to be resolved well enough. Then, there 
is the large length scale of the supermassive black hole, 
which sets the distance to the wavezone, where the grav- 
itational wave signal is to be extracted. Finally, since a 
typical EMRI source for LISA will make about 10 4 - 10 5 
orbits, long-term evolutions will be required to produce 
gravitational wave templates of use to data analysis. 
These considerations conspire to make EMRI modeling 
a difficult problem for numerical relativity. 



*URL: http:/ /www. phys.ufl.edu/ift/ 
tURL: htt p://www.cct.lsu.edu/! 
*URL: http://relativity.phys.lsu.edu/ 



At some point, the steady advance of computational 
technology will allow the numerical relativity commu- 
nity to tackle the full dynamics of this binary system. In 
the meantime it seems prudent to develop approximate 
schemes that will reliably produce templates of adequate 
accuracy. One such scheme takes advantage of the small 
mass-ratio (say ii/M) and treats the problem in a pertur- 
bative fashion. At lowest order, the system is but an in- 
finitesimal test mass moving in a black hole spacetime, for 
which we know the motion to be geodesic in that space- 
time. For such the mature formalism of black 
hole perturbation theory is able to accurately calculate 
first-order metric perturbations, from which one infers 
gravitational waveforms. The foundations for these sorts 
of calculations [1, 0, Q were laid out beginning over fifty 
years ago. However, the accuracy requirements of LISA, 
particularly on the phase of the waveform throughout 
the entire mission lifetime, demand that our models go 
beyond this leading order case. (An instructive scaling 
argument can be found in Sec. 11.1 of [(|). One thus has 
to consistently take into account next-to-leading order 
effects on the motion of the particle and the waveform. 
From this perspective, the dynamics of the inspiral are 
viewed as the motion of a finite (but small) point mass 
in a background black hole spacetime. The goal is then 
to determine these finite-mass effects on both the motion 
of the point mass and the corresponding gravitational 
waveform. 

These effects are attributed to the self-force; the dissi- 
pative part of which is the more familiar phenomenon of 
radiation reaction. For a point mass moving in flat space- 
time, this effect is entirely local. Only emission at a given 
instant affects the motion of the particle at that same in- 
stant. In curved spacetime, due to scattering with the 
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curvature, the motion of the particle is affected by fields 
it gave rise to in its causal past. This makes the resulting 
physics much richer than of the flat spacctime case. The 
same scenario also applies to scalar and electric charged 
particles moving in curved spacetime. While often inter- 
esting in their own right, these also serve as useful and 
technically less-demanding toy problems for testing new 
methods and techniques. 

The effect of the self-force is a small acceleration on 
the point mass resulting in a secular deviation away from 
what would otherwise have been geodesic motion in the 
background spacetime. This self-force will need to be cal- 
culated and used to modify the motion of the point mass 
as often as is practical throughout the course of the inspi- 
ral. Formal expressions for this self-force given in terms 
of an integral over the particle's entire past history have 
been ironed out in the literature @, @, ls| , but these are 
hardly convenient for practical calculations. The chal- 
lenge is then to come up with an efficient way to compute 
self-forces based on these formal expressions. Several use- 
ful prescriptions have been developed to address this is- 
sue fiolll l|. but most influential of these is the mode-sum 
prescription [f2l . [l3l . [l4| , which we shall describe below. 

An alternative viewpoint of the self-force scenario is 
that instead of a backreacting "force" accelerating the 
point mass in the background spacetime, the motion of 
the point mass is really geodesic motion on the distorted 
background geometry [g, [l5| . In this framework, the 
task is to determine the appropriate distorted geometry 
upon which to impose geodesic motion (or equivalcntly, 
the correct smooth potentials governing the motion of 
scalar and electric charges) . The procedure thus involves 
first computing the metric perturbation h a b induced by 
the point mass on the background spacetime metric g a b 
(which would be divergent at the location of the point 
mass), and then appropriately regularizing this metric 
perturbation to give the correct smooth perturbation h^ b . 
The motion of the point mass will then be geodesic mo- 
tion in the perturbed spacetime g a b + h^ b . 

This perspective has been useful on both theoretical 
and practical fronts. Most notably, it has reconciled the 
notion of a self-force with our understanding of the equiv- 
alence principle [l^|. It has also served as the basis of 
convenient variants of the ori gina l mode-sum prescrip- 
tion [U, EE EE EI HI M, HI M Hi , and has thus ad- 
vanced our understanding of this important calculational 
technique. The method presented in this manuscript is 
another off-shoot of this alternative viewpoint. 

Much progress has been made in the calculation of the 
self-force on a charge that moves momentarily along some 
prescribed geodesic of the background spacetime. In par- 
ticular, the mode-sum paradigm has contributed tremen- 
dously to our understanding of the elements of a self- force 
calculation and continues to serve as the conceptual back- 
drop upon which all other calculational schemes are to 
be understood. It has now been employed to evaluate 
the self-force or self-force effects in a variety of contexts - 
scalar [H, H3, d, [H, HH, [H, HI] , electromagnetic [H, 



and gravitational 0, [2?], HI] - for point sources moving 
along geodesies in a Schwarzschild background. Results 
for the Kerr spacetime are rare. The first calculation of 
a self force on a scalar charge moving through this back- 
ground has been achieved only recently [1^] . Despite this 
good progress, however, little has been done to achieve a 
dynamic calculation of the self-force that is then used to 
implement backreaction on moving point sources. 

One of the reasons for this is that computing a self- 
force is a complicated process. A typical mode-sum cal- 
culation first requires a decomposition of the problem 
(i.e. fields and sources) into modes with, say, spherical 
or spheroidal harmonics. This is done in order to avoid 
having to handle the divergence in the physical retarded 
field numerically. Each mode component of the retarded 
field turns out to be finite at the location of the charge, 
and thus, numerically accessible. It is the sum of these 
modes that diverges. Buried in each mode is the piece 
that actually contributes to the self-force. The central in- 
sight of the mode-sum prescription is a way to access this 
relevant piece, based on an asymptotic analysis of the di- 
vergence in the retarded field. In calculating the self-force 
then, each computed mode is appropriately regularized 
(using an analytic expression determined by the asymp- 
totic analysis) to leave the piece that contributes to the 
self-force. The regularized pieces are then summed to get 
the full self-force. Convergence of this sum is typically 
slow, going as ~ f/^™, where I is the maximum mode 
number, and n is determined by the degree to which one 
analytically characterizes the asymptotic behavior of the 
divergent physical fields. 

In [23| (henceforth referred to as Paper I), we intro- 
duced a method for calculating the self-force designed 
principally to obviate the apparent necessity of a de- 
composition into modes in order to perform the regu- 
larization of the retarded field. Through this alterna- 
tive method, we proposed a (3+1) approach to self-force 
calculations, which fits squarely with the expertise and 
infrastructure found in the numerical relativity commu- 
nity. Our technique is reviewed in better detail below. 
Its core idea, however, is straightforward: rather than 
regularizing the retarded field, one can instead appro- 
priately regularize the source term in the field equation 
from the start, and thus have the evolution codes deal 
with sufficiently-differentiable fields and sources that re- 
quire no further rcgularization. In other words, we deal 
with regularization not as a post-processing step, but as 
preliminary work that needs to be performed before any 
numerical run. This work consists of appropriately re- 
placing the delta-function representation of a point mass 
source by a regularized effective source. Designing this 
effective source for the wave equation can be done in 
such a way that a numerical evolution yields a differen- 
tiable field whose gradient at the location of the particle 
automatically gives the full self-force. In addition, this 
resulting regular field is such that it becomes the physical 
retarded solution in the wavezone, from which fluxes and 
the all-important waveforms can be extracted. 
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Two important features of this approach stand out. 
First, the evolution code never has to deal with diver- 
gent quantities (though the fields and sources are of finite 
differentiability); and second, both self-force calculation 
and waveform extraction are trivial, with accuracies lim- 
ited chiefly by the accuracy that can be provided by the 
evolution code. It is these two features that call out to the 
numerical relativity community at large, for in effect all 
that is required for a self-force calculation is a (3+1) code 
that can accurately evolve wave equations with sources of 
limited differentiability. The ease with which one calcu- 
lates a self-force in this approach suggests no significant 
impediment to implementing backreaction on the parti- 
cle. 

The idea of finding a good substitute for the delta 
function source in the context of self-force calculations is 
also being pursued by others in similar ways [13, [U, HI] . 
Barack and Golbourn [30l | introduced a technique consist- 
ing of an m-mode (azimuthal mode) regularization of the 
delta function source. Their regularization is also guided 
by analysis of the singular behavior of the retarded field 
at the location of the particle. First one solves (2+1) 
wave equations with regularized sources, then extracts 
the contribution to the self-force due to each azimuthal 
mode, and finally sums these to get the full self-force. 
Their approach is similar to ours in that it provides a 
way of representing point sources on a grid, but in keep- 
ing with the general strategy of the original mode-sum 
procedure (which was an i-mode sum), it is likely to in- 
herit some of the properties our approach seeks to over- 
come. 

A new approach to evaluating the retarded field by 
Canizares and Sopuerta [32l | splits the problem into in- 
ner and outer domains marked by the location of the 
point charge. They situate the point source along the 
boundary shared by both domains, and then just impose 
appropriate jump conditions on the fields that cross the 
boundary. With the benefit of having to deal with only 
smooth fields, this method takes advantage of the expo- 
nential convergence of a pseudospectral implementation 
for evaluating the retarded field. In calculating a self- 
force, however, they are still restricted to performing a 
mode-sum of what remains after regularizing the output 
of their evolution code. A chief advantage of our method 
is precisely the fact that it escapes the requirement of a 
mode decomposition. 

In Paper I, we reported an implementation of our 
method using a time-domain (1+1) code to compute the 
self-force and retarded field in the wavezone. By first 
breaking into modes, this implementation ran counter 
to the very motivations underlying our technique. This 
was done, however, mainly to provide a quick proof-of- 
principle, and to establish a more direct connection with 
more familiar approaches to self-force calculation. 

In this work, we report for the first time on the fea- 
sibility of our technique in its intended setting. As a 
result, we have achieved the first calculation of a self- 
force in a (3+1) framework. Two different codes [33l. l3~j| 



were employed in the implementation of our method, 
both originally intended for numerical relativity appli- 
cations. As in Paper I, we focus on the simplest possible 
strong-field scenario involving a scalar point charge inter- 
acting with its own scalar field while in a circular orbit 
around a Schwarzschild black hole. In choosing to deal 
with this simple case, we also intend for this document 
to be a self-force primer and an invitation to numeri- 
cal relativists who are on the lookout for new challenges. 
The insights gained from self-force analyses should prove 
useful to those wishing to tackle the extreme-mass-ratio 
regime of black hole binaries. 



Outline and notation 

Section [II] describes our formalism for replacing a 
point-particle delta-function source in a standard field 
equation with a particular abstract effective source. The 
effect is that not only does the resulting field equal the 
usual retarded field in the wave zone, but also the field 
is finite and diffcrentiablc at the particle. This allows it 
to be used directly in calculating the self-force acting on 
the particle. Some details regarding the effective source 
are given in Sect. IIIIl 

We describe a practical test application for our ap- 
proach to self- force computations in Sect. HVl This 
test application has been previously well studied by the 
self-force community by using more traditional self-force 
techniques which are not particularly adaptable to nu- 
merical relativity [H, [H IS El • 

Section [V] contains the implementation details of mod- 
ifications, with an eye on applications to self- force prob- 
lems, of two different previously developed numerical rel- 
ativity projects. 

Details of the results following from the applications of 
these two codes to our test are given in Sect. I VII We com- 
pare the time and the radial components of the self-force, 
calculated by our two independent codes, with each other 
and with very accurate (but tediously obtained) well- 
known frequency-domain results [l|| Qj], EH, El- The 
time component of the self-force removes energy from 
the particle, and we also check its consistency with the 
energy flux via radiation down the black hole and out at 
infinity. 

Section [VI 1 1 gives a summary of the apparent strengths 
and weaknesses of our effective source method for regu- 
larizing self- force problems. In very general terms we 
describe how currently available computer codes might 
be adapted specifically to self-force problems. 

We have three appendices. Appendix [Al gives a 1+1D 
example of applying traditional finite differencing oper- 
ators to a wave equation where the source is of limited 
differentiability. This elucidates the discussion of con- 
vergence. In Appendix [B] we derive the relationship be- 
tween between the time component of the self-force and 
the radiative energy-flux into the black hole and out at 
infinity, in the case of a circular orbit and a scalar field. 
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And wc describe a very elementary, illustrative flat-space 
toy problem in Appendix [U] which demonstrates how a 
problem involving a delta-function point source can be 
transformed into one with a smooth source in a mathe- 
matically precise way. 

For our tensor notation, we denote regular four- 
dimensional space time-time indices with letters taken 
from the first third of the alphabet a,b, . . . ,h , indices 
which are purely spatial in character are taken from the 
middle third, i,j,...,q and indices from the last third 
r, s, . . . , z and also 9 and <p are associated with particular 
coordinate components. The operator V a is the covariant 
derivative operator compatible with the metric at hand. 
Partial derivatives with respect to t are denoted dt, and 
with respect to a generic spatial coordinate by cV 

For the Schwarzschild metric, we use a coordinate sys- 
tem introduced by Eddington [35| and commonly known 
as Kerr-Schild coordinates to describe a Schwarzschild 
black hole. 

Our use of the 3+1 formalism follows York [36| in all 
aspects except his labels for tensor indices. 



II. FIELD REGULARIZATION FOR A SCALAR 
CHARGE 

In this section we review the discussion of our method 
found in Paper I. We shall discuss it for the case of 
a scalar point charge moving in curved spacetime. A 
typical self-force computation first involves solving the 
minimally-coupled scalar wave equation with a point 
charge q source, 

^ = -4*9 I S -^M dr, (1) 

for the retarded field ?/i ret . Here V a is the derivative 
operator associated with the metric g a b of the background 
spacetime and 7 is the world line of the charge defined 
by z a (r) and parameterized by the proper time r. The 
physical solution of the resulting wave equation will be a 
retarded field that is singular at the location of the point 
charge. A formal expression for the self-force given by 

F a (r) = q(g a b + u a u b )V b ^ c \z{r)) (2) 

would thus be undefined without a proper regularization 
prescription. Early analyses @, H, Ht| were based upon a 
Hadamard expansion of the Green function, and showed 
that for a particle moving along a geodesic the self-force 
could be described in terms of the particle interacting 
only with the "tail" part of ip let , which is finite at the 
particle itself. The mode-sum prescription is effectively 
a way of regularizing the righthand side of Eq. ^ to re- 
trieve the force due to this tail part. Later it was 
realized that a singular part of the field ^ s which exerts 
no force on the particle itself could be identified as an 
actual solution to Eq. ([TJ in a neighborhood of the par- 
ticle. A formal description of ip in terms of parts of the 



retarded Green's function is possible, but generally there 
is 110 exact functional description for -0 s in a neighbor- 
hood of the particle. Fortunately, as is shown in [l(|, 
an intuitively satisfying description for ijj results from a 
careful expansion about the location of the particle: 

^> s = q/p + 0{p z lK A ) as p -> 0, (3) 

where 1Z is a constant length scale of the background 
geometry and p is a scalar field which simply satisfies 
p 2 = x 2 + y 2 + z 2 in a very special Minkowski-like locally 
inertial coordinate system centered on the particle, first 
described by Thorne, Hartle and Zha ng [38l. l39j and ap- 
plied to self- force problems in Refs. [til. [ml4Cj] . A detailed 
discussion of these coordinates can be found in [|| [l6| 
and Appendix A of Paper I. Not surprisingly the singu- 
lar part of the field, which exerts no force on the particle 
itself, appears as approximately the Coulomb potential 
to a local observer moving with the particle. 

Our proposal for solving Eq. fTj), and determining the 
self-force acting back 011 the particle now appears ele- 
mentary. First we define 

V? = q/p (4) 

as a specific approximation to ip s . By construction, wc 
know that ^ S is singular at the particle and is C°° else- 
where. Also, within a neighborhood of the world line of 
the particle 

V a V a ^ s = -47T9 / ^%^rfr + 0(p/K 4 ), 
J-t V~9 

as p — > 0. (5) 

It must be pointed out that for Eqs. and (O to be 
valid, the Thorne-Hartle-Zhang (THZ) coordinates must 
be known correctly to 0(p 4 /7Z 3 ). Knowing the THZ co- 
ordinates only to 0(p 3 /lZ 2 ) would spoil the remainder 
in Eq. ([5]) which would then have a direction dependent 
discontinuity in the limit as p — » 0. The local coordinate 
frame must be known precisely enough in terms of the 
global coordinates for the Coulomb-like potential to be a 
good representation of the local singular field. 

Next, we introduce a window function W which is a 
C°° scalar field with 

W = 1 + 0{p 4 /K 4 ) as p -> 0, (6) 

and W — > sufficiently far from the particle, in particu- 
lar in the wavczone and at the black hole horizon. The 
requirement that 11" approaches 1 this way, i.e. 0(p l ). is 
explained below. 

Finally wc define a regular remainder field 

V^ R = V/ Ct - w^ s (7) 

which is a solution of 

v°VaV' R = -v a v Q (py^ s ) 

-47rg / - v " dr (8 
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from Eq. (JTJ) - Because of our use of ip s as an approxi- 
mation of the full singular field, ip R will be contaminated 
with 0(p 3 /7?. 4 )-pieces that are only C 2 at the location of 
the charge but do not affect the self-force. 
The effective source of this equation 



S eff =-V Q V a (W^ s )-4^ 



-9 



dr (9) 



is straightforward to evaluate analytically, and the two 
terms on the right hand side have delta-function pieces 
that precisely cancel at the location of the charge, leaving 
a source which behaves as 



S eS = 0(p/n A ) as p -> 0. 



(10) 



Thus the effective source S c g is continuous but not neces- 
sarily differentiable, C°, at the particle while being C°° 
elsewhere 1 . 

A solution ?/> R of 

V a V Q ^ R = S cS (11) 
is necessarily C 2 at the particle. Its derivative 
V a i/; R = V a (# ct - W^ s ) - i> s V a W 

= v a (^ t -^ s ) + o(p 2 /n i ) P ^o (12) 

provides the approximate self-force acting on the particle 
when evaluated at the location of the charge. It should 
be clear why the behavior of W is chosen as in Eq. © : A 
window function with this behavior would not spoil the 
0(p 2 /H )-error already incurred by using the q/p ap- 
proximation for the singular field. Also, in the wavezone 
W effectively vanishes and ip R is then identically -0 rct and 
provides both the waveform as well as any desired flux 
measured at a large distance. 

General covariance dictates that the behavior of S e s in 
Eq. @ may be analyzed in any coordinate system. But, 
only in the specific coordinates of Refs. [38[ and (3j|, or 
the THZ coordinates, is it so easily shown [l6| that the 
simple expression for ip° in Eq. ([3]) leads to the 0(p/lZ 4 ) 
behavior in Eq. (|10|) and then to the C 2 nature of the 
solution ip R of Eq. (JTTJ) . 

Self-consistent dynamics of a scalar charge requires 
that the self-force act instantaneously. Thus a simul- 
taneous solution of the coupled equations 

V a V a V R = S cS (x(t),u(t)) (13) 
dv b 

= q(g bc +A C )V/ (14) 

dr 

evolves ip R while self-consistently moving the charge via 
the self- force of Eq. fl4|). 



Our method effectively regularizes the field itself rather 
than the gradient of the field, and this regularization is 
implicit in the construction of the effective source, as 
opposed to most existing self-force calculations in which 
the divergent pieces of the individual modes are explicitly 
subtracted out. Once ip^ is determined, our method has 
no need for any further regularization. The derivatives 
of V> R determine the self-force, providing instantaneous 
access; while ip K is identical to ip TQt in the wavezone, 
allowing direct access to fluxes and waveforms. 

The tedious aspects of our method reside primarily in 
the construction of the effective source. This is mainly 
due to the need for the transformation from THZ coor- 
dinates to the background coordinates. This transfor- 
mation is a function of the location and four-velocity of 
the particle at any given instant. For fully consistent 
dynamics where the particle location and four-velocity 
are constantly being modified, this transformation will 
itself be changing. Thus, this coordinate transformation 
will unavoidably have to be determined and then applied 
numerically. 

Once the effective source is appropriately constructed, 
the only remaining requirement is a code capable of 
evolving the wave equation with a C° source. 



III. EFFECTIVE SOURCE 

At the heart of our approach is the use of a convenient 
regular representation of a point particle source. We refer 
to this as the effective source. The two main elements 
which enter this are (1) the approximate singular field, 
■0 s = q/p, whose explicit form in terms of the chosen 
background coordinates depend on the position and four- 
velocity of the particle, and (2) the window function, 
W , whose main purpose is to localize the support of the 
approximate singular field to within the vicinity of the 
particle. 

In tackling the same physical test application as in Pa- 
per I, no modifications of ip were needed for our (3+1) 
runs, apart from a trivial replacement of the background 
coordinates in which to express the effective source. 
(Here we use ingoing Kerr-Schild coordinates as opposed 
to the Schwarzschild coordinates of Paper I) . 

However, for the current implementation, we did seek 
out a more adaptable window function. In (1+1), it 
proved sufficient to use a simple window function hav- 
ing a Gaussian-like profile in r: 



W(r) = cxp 



(r - r Q ) 



N 



(15) 



1 With p 2 = x 2 + y 2 + z 2 , a function which is 0(p n ) as p — > 0, is 
at least C n where p = 0. 



where r is the radius of the circular orbit in Schwarz- 
schild coordinates, while N and a are parameters to be 
chosen according to the requirements described in SjH] It 
is easily verified that all of these required conditions can 
be met for a sufficiently large N. In principle, these con- 
ditions make it a reasonable choice regardless of the nu- 
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merical implementation. In practice, however, this orig- 
inal effective source had some properties that could po- 
tentially burden certain (3+1) codes. (Some results from 
our early runs with the original source were, in fact, what 
motivated the construction of a new one.) Specifically, 
the choice of a Gaussian-like window leads to significant 
large-amplitude, short-scale (~ a) structure away from 
the particle. This was not an issue in the (1+1) case, 
where high /'-resolution (Ar ~ Af/25) and high angular 
resolution (with a spherical harmonic decomposition go- 
ing to as high as L = 39) were practical. Of course, the 
extra structure away from the particle need not necessar- 
ily be a problem for all (3+1) codes. One can maintain 
the Gaussian-like window and simply adjust its width a 
to lessen the artificial short-scale structure. Some of the 



runs presented below were performed with this original 
window, using N = 8, r a = 10M and a = 5.5M. The 
width was chosen in order to make the profile as wide as 
possible while still effectively vanishing before the hori- 
zon is reached. These runs show sufficiently good results 
as well. 

Nevertheless, there is merit in using a more flexible 
window function; for instance, one with more adjustable 
parameters that can be tuned to the needs of any (3+1) 
code. A convenient choice makes use of the smooth tran- 
sition functions introduced in |4lj |. Like in Paper I, we 
have chosen to apply a window function only along the 
r-direction, in keeping with the spherical symmetry of 
the background spacctime. 



J 



Consider the smooth transition function 

r o, 



f(x\x ,w,q,s) 



1 1 ,1 s , 

= < — I — tanh — < tan 

2 2 7T 



1 1. 



2^(z-Zo) 



tan 



2w 



(x - x ) 




This is a function that smoothly transits from zero to 
one in the region xq < x < xq + w. It comes with four 
adjustable parameters {xq, w, q, s}: 

1. Xq: defines where the transition begins. 

2. w: gives the width of the transition region. 

3. q: determines the point = xo + (2w/ir) arctang 
where the transition function f(x) = 1/2 . 

4. s: influences the slope s(l + q 2 )/(2w) at Xi/ 2 after 
w and q are chosen. 

Using this transition function, a window function for a 
particle at r = R could be 

Vf/M - J f( r \( R - S l -wi),wi,qi,si) r<R , 
W[ ' ~ \ l-f(r\(R + 5 2 ),w 2 ,q 2 ,s 2 ) r>R 

and W(r) = 1 in the region R — Si < r < R + S 2 . 

This satisfies all of the key requirements for a window 
function (and more): 

(a) W(R) = 1; 

(b) d n W/dr n \ r = R = 0, for all ra; 

(c) W = if r G [0, R - Si - wi] U[R + 5 2 + w 2 ,oo) 
(thus making it truly of compact support); 

(d) and W = 1 if r G [{R - Si), {R + S 2 )\- 



For the actual runs that used this window function, we 
settled on the following choices for these parameters: 
{Si = S 2 = 0M;<7i = 0.6, 92 = 1.2; si = 3.6, s 2 = 
1.9; wi = 7.9M, w 2 = 20M}. The inner width wi was 
chosen so that the window and effective source go to ex- 
actly zero just outside the event horizon. The rest were 
picked after extensively looking at many parameter com- 
binations. The primary criteria were simply that the 
effective source would be sufficiently small everywhere 
and that it did not possess structure at extremely small 
scales. A systematic search for the optimal set of param- 
eters vis-a-vis its effect on self-force accuracy was not 
conducted in this study, and is left for future work. 

One important attribute of the new window function 
is that, for a wide range of parameter choices, it leads 
to an effective source whose over-all structure away from 
the particle is significantly less pronounced than that pro- 
duced by the original Gaussian-like window. Comparing 
the new effective source in Fig. [T] with the original source 
used in Paper I (shown as Fig. 1 of that paper), one 
notes immediately that the artificial structure resulting 
from the new window is almost two orders of magnitude 
smaller. Moreover, this structure is mainly located at 
r <R (where R = 10M). 

It is instructive to look at the structure of S e fi at the 
location of the particle. The effective source, Soft, is C° 
at the particle due to the level of the approximation used 
for the singular field tp s . This C° behavior is sufficient 
for calculating the self-force. In our approach, this yields 
an evolved regular field -0 R that is C 2 at the location of 
the charge, from which derivatives can be computed to 
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Figure 1: Equatorial profile of the new effective source, S e ff, 
at t = 0. The axes are defined simply by x — rsin9cos<f) 
and y = r sin 6 sin cj>, where r, 9, <j) are just the Schwarzschild 
coordinates (or the Kerr-Schild coordinates of Sec. V A). The 
charge in this plot is located at X = 10M and Y = 0, where 
the C° behavior of the source is not apparent on this scale. 
Note that much of the structure induced by the new window 
function is between the charge and the event horizon. 



give the self-force. In Fig.[2]the C° nature of the effective 
source is revealed. The effective source is certainly a non- 
singular representation of a point charge source which is 
amenable to the (3+1) codes we have used for calculating 
the self-force. 




Figure 2: Sgg zoomed in at the location of the charge. 



IV. TEST APPLICATION 

The physical scenario that is analyzed in this paper 
involves a particle with mass m and scalar charge q in 
a perpetual circular orbit around a Schwarzschild black 



hole while emitting scalar radiation. The effects of the 
self-force (which include the effects of the emission of 
radiation) would lead to the gradual decay of the circular 
orbit. For simplicity in this analysis we keep the charge 
in a circular orbit and compute the external force needed 
to counteract the scalar self- force. 

With the charge in perpetual circular motion, and in 
the absence of other external sources which may violate 
this symmetry, the system is helically symmetric. For 
any field G, there must then exist a helical Killing vector 
£ a , such that 

£ C G = 0. (18) 
In Schwarzschild coordinates, this Killing vector is simply 



d 



d_ 

at 



o 



d_ 



and 



(19) 



(20) 



For the circular orbit problem we have chosen, the four- 
velocity u a is tangent to the Killing field at the location 
of the particle. Thus, V a <^> is already orthogonal to u a , 
and the self-force is given directly by 



(21) 



with no need for the projection operator, present in 
Eq. (|14[) . for our test application. There arc then only 
two independent components of the self-force, F t and F r , 
with Fj, = -F t /Q from Eq. ([20]). and Fg = 0, by virtue of 
the system being reflection symmetric about the equator. 

For circular orbits, there exists a useful relation be- 
tween the scalar energy flux and F t . In terms of the 
Kerr-Schild coordinates described next, this appears as 



dE 
~dt 



where 



dE 
~dt 



dE 
~dt 



=2JVf 



r=R 



= -,/!- ™F., 



4M 2 <* i? dfl, 



(22) 



(23) 



and 



dE 
~dt 



R' 



r=R„ 



R 



R-2M 



2M 

IT* 



2AI 

— ) i\,d r i\) 



dQ. 



(24) 



Here, r is the radius of the circular orbit, and R is the 
finite outer extraction radius. The field ip is actually the 
retarded field, but one can instead use as long as 
the surface integrals are evaluated outside the support of 
the window function, where (by design) ■0 ret = This 
simple relation is proved explicitly in Appendix [B] We 
use this as a consistency check on our self-force results. 
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A. Coordinates for a Schwarzschild black hole 

We describe the Schwarzschild metric as 

9ab = Vab + Hk a k b (25) 

using a coordinate system first identified by Eddington 2 
and commonly known as Kerr-Schild ingoing coordinates 
(t, x, y, z), where r) a b is the flat Minkowskii metric with 
(—1, 1, 1, 1) along the diagonal, 



k a dx a = —dt 
= -dt 



x , y , z i 

— dx-\ — dy H — dz 
r r r 



dr 



which is the ingoing principle null vector, and 



(26) 



(27) 



with r 2 — x 2 + y 2 + z 2 . This is equivalent to the usual 
Schwarzschild form of the metric 



ds 2 



7M\ 



dt 1 



dr 2 



r J \-2Mjr 
+ r 2 (d9 2 + sin 2 6 dcj) 2 ). (28) 



with the Schwarzschild time coordinate t related to the 
Kerr-Schild coordinates by 



i = t-2Mln(r/2M - 1) 



(29) 



and the usual flat space relationships between (r, 9, <p) 
and (x,y,z). 

The Kerr-Schild form of the metric is popular in the 
numerical relativity community because a constant t hy- 
pcrsurface is non-singular and horizon penetrating which 
allows for convenient imposition of boundary conditions 
or for excision. 

However, it can be confusing to compare the com- 
ponents of the self-force as evaluated in these Kerr- 
Schild coordinates with the components as evaluated in 
Schwarzschild coordinates. The Kerr-Schild radial co- 
ordinate tks equals the Schwarzschild radial coordinate 
^Sch, but constant-i and constant-i surfaces are not the 



2 Actually Eddington [35l | and Finkelstein [42l | wrote down the 
Schwarzschild metric using the precise coordinates of Eq. 1261 . 
But, somehow the Eddington-Finkclstein duo arc associated with 
a coordinate system that contains an ingoing or outgoing null co- 
ordinate, although neither explicitly introduced or used such a 
null coordinate. While Kerr and Schild (nearly forty years after 
Eddington) described the Kerr metric in a form that reduces to 
Eq. H26j l in the Schwarzschild a — » limit. Bowing to current 
conventions of the numerical relativity community rather than 
to historical accuracy, we label the coordinate system in use as 
"Kerr-Schild." 



same. The relationships between the components of the 
self-force for these two coordinates systems are 



^Sch _ ^KS 



^Sch _ 



2M 



r n - 2M 



(30) 
(31) 



B. The 3 + 1 version of the Schwarzschild metric in 
Kerr-Schild coordinates 

For the Schwarzschild metric in Kerr-Schild coordi- 
nates the contravariant form of the metric (1231) is 



g ab = r] ab - Hk a k b 



(32) 



With the 3 + 1 formalism [36j the contravariant compo- 
nents of the metric are closely related to the lapse func- 
tion a, shift vector (3 l and spatial metric 7 y of a foliation 
of spacctime by 



-a~ 2 PI a 2 

J3 l /a 2 - f3 l ft/a 2 

This relationship gives 

- g u = 1 + H = a' 2 , 

g lt = Hx l /r = (3 l /a 2 , 

and 
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which implies that 



7 u = rfj 



H x l x J 



(33) 

(34) 
(35) 

(36) 
(37) 



l + H r 2 

Also the determinants of the metrics arc related by 

V^g = "V7- (38) 

C. The wave equation in Kerr-Schild coordinates 

A form of the wave operator convenient for computa- 
tion is 

1 



-9 



--d a 



99 ab d b ^ 



(39) 



In the 3 + 1 formalism, after substitution for the con- 
travariant components of the metric from Eq. (|32p and 
with the time- independence of g ab , we obtain 

a 2 V Q V a ^ = -cW + ^cW 



V7 \« 



a 



V7 



aV7( 



(40) 



for the wave equation in the Schwarzschild geometry 
with Kerr-Schild coordinates and a, (3 l and y lJ given in 
Eqs. (EU-iU). 



9 



V. (3+1) IMPLEMENTATIONS 

In the following we describe the finite differencing and 
pscudospectral codes used in the numerical experiments. 



3D multi-block finite difference code 



We solve the wave equation ([13]) for tp on a fixed 
Schwarzschild background with a source over a multi- 
block domain using high order finite differencing. The 
code is described in more detail in [33[, here we will just 
summarize its properties. We use touching blocks, where 
the finite differencing operators on each block satisfies a 
Summation By Parts (SBP) property and where charac- 
teristic information is passed across the block boundaries 
using penalty boundary conditions. Both the SBP opera- 
tors and the penalty boundary conditions are described in 
more detail in [4311 . T he code has been extensively tested 
and was used in [44| to perform simulations of a scalar 
field interacting with Kerr black holes and was used to 
extract very accurate quasinormal mode frequencies. 

After the standard 3+1 split, the wave equation is 
written in first-order in time, first-order in space form 
in terms of the variables p = dtip and fa = diip. The 
system of equations being integrated is then 



dtp 
d t fa 



y^f L V Oi 1 



a 2 S cS} 



dip, 
P- 



(41) 
(42) 
(43) 



where Eq. fllT]) follows from Eq. (J40J) with g ij = j ij - 
a~' 1 (5 % ^ as in Eq. (|3l)|) . and Eq. (|42[) is an elementary 
consequence of the definition of (f>i. The primary dy- 
namical variables are u = (p,fa), while ip is evolved via 
an ordinary differential equation (no spatial derivatives). 
Across a boundary with unit normal vector £j the char- 
acteristic modes are: 



(44) 
(45) 



fa - ¥(f>j£i 



Where the speeds of the two transverse modes in Eq. 
are A = and the speeds of the two normal modes in 
Eq. gS]) are A ± = ± a. 

The only necessary modifications to the code described 
in [33] , in order to apply it to the problem at hand, were 
the addition of the source term in Eq. (|4"Tj) and to add 
code to interpolate the time derivative p and the spatial 
derivatives fa of the scalar field to the location of the 
particle. 

In addition some optimizations were performed. 
OpenMP pragmas and directives were added to allow for 
simultaneous OpenMP and MPI parallelization for better 
performance on modern multi-core machines. Also a load 



balancing issue arose that could potentially lead to very 
poor scaling because S c s is expensive to calculate only in 
the spherical shell where it is non-zero. This issue was 
solved by adding data structures that were distributed 
evenly among all MPI processes, with just the right size 
and shape to cover the spherical shell. The source is then 
evaluated first (all processors working simultaneously) on 
this distributed data structure and then copied into the 
main 3D grid functions. 



1. Boundary conditions and initial data 

The simulations below were all performed using the 6- 
block system, providing a spherical outer boundary and 
spherical inner excision boundary without any coordi- 
nate singularities. We use the Schwarzschild solution 
in Kerr-Schild coordinates as the background metric for 
the scalar field evolution. The inner radius was chosen 
to be i?i n = 1.8M and the outer boundary was chosen 
to be at i?out = 400M in most cases (it was placed at 
Rout = 600M in a few runs for more accurate extraction 
of the fluxes). Since we are using SBP finite differenc- 
ing operators we can evaluate the right hand sides for 
the evolution equations, i.e. dtu = (dtp, dtfa), even as we 
approach the outer boundary (using more and more one- 
sided stencils). At the outer boundary we then convert 
both u and dtu to characteristic variables using Eqs. ([44]) 
and 1(15]). i.e. we obtain {w^w^ and (d t w^,d t w ± ). We 
only have to apply a boundary condition to w~ , since 
this is the only incoming mode. We do this by adding 
a suitable penalty term (only at the outer boundary) to 
dtw~ of the form T(g — w~), where T is a penalty pa- 
rameter that has to be chosen to be consistent with the 
SBP operator and the speed of the mode in order to 
achieve stability (see more details in (33|) and g is the 
desired incoming characteristic mode. In this case we 
use g = 0, i.e. zero incoming mode. We then transform 
{dtW®, d t w + , d t w~ +T(g — w~)) back to a new d t u that is 
used by the time integrator to update the primary vari- 
ables. At the inner boundary, the geometry ensures that 
all characteristics leave the computational domain; i.e. 
there are no incoming modes and therefore we do not 
apply any boundary condition there. 

We do not, a priori, know the correct field configu- 
ration and start the simulation with zero scalar field 
ip(t = 0) = 0, zero time derivative pit = 0) = and 
zero spatial derivatives fa(t = 0) = 0, as if the scalar 
charge suddenly materializes at t — 0. After the system 
is evolved for a few orbits, the initial transient has de- 
cayed and the system approaches a helically symmetric 
end state. We used the 8-4 diagonal norm SBP oper- 
ators and added some compatible explicit Kreiss-Oligcr 
dissipation to all evolved variables. 

With this code we have performed runs for a scalar 
charge on circular orbits of radius r Q = 10T\7 with 
both the wide Gaussian profile window (TV = 8 and 
a = 5.5M) and the smooth transition function win- 
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dow (Si = 8 2 = 0M;qi = 0.6, q 2 = 1.2; Si = 3.6, s 2 = 
1.9; wi = 7.9M,w 2 = 20M). We find that the extracted 
self-force is independent of the window function (as it 
should be) and that the only difference between the runs 
is in the shape and amplitude of the initial scalar wave 
pulse. 



2. Convergence 

The convergence of the code has been extensively 
tested in [43| where the evolution of a plane wave moving 
across a spherical grid was used as a test problem (i.e. no 
source). It was shown that for all implemented finite dif- 
ferencing order the code was converging at the expected 
order. For example for the 8-4 SBP operators used here 
we found the expected fifth order global convergence. 

As the source is only C° at the world line of the par- 
ticle, it is to be expected that the scalar field will be C 2 
there while the main evolution variables p = dip/dt and 
4>i = dtp/dx* should be C 1 on the world line of the parti- 
cle and C°° everywhere else. With the finite differencing 
code, there will be some stencils which are penetrated 
by the worldline at a particular event. For those sten- 
cils, the finite differencing errors will be affected by the 
limited differentiability of both the source and the field 
at the particle. We would expect that any traditional 
centered finite differencing operator applied to a C 1 field 
(regardless of order) should then only be first order ac- 
curate: the second derivative is discontinuous across the 
world line and so the second order terms in the Taylor 
expansion of the operator will not cancel. 

Naively one would then expect that the solution for p 
and cf>i at the particle and thus the extracted self-force 
would only converge to first order. However, as shown in 
Appendix [X] for the wave equation in 1+1D, the errors in 
p in fact converge at second order in the L2-norm for a C° 
source. In the Appendix, it is also shown that the error is 
of high frequency with the frequency increasing with res- 
olution. Thus, for our test application we cannot demon- 
strate pointwise convergence for the quantities p and <fii . 
But we expect that the amplitude of any noise generated 
near the particle location will converge at second order. 
We find below that the extracted self-force components 
at the location of the particle are indeed noisy, but that 
the noise converges to zero at second order. 



B. The pseudospectral code 



We solve the wave equation (|T3[) for tp on a fixed 
Schwarzschild background with a source using pseu- 
dospectral techniques. 

We use the SGRID code |U HI to numerically 
evolve V- This code uses a pseudospectral method in 
which all evolved fields are represented by their values 
at certain collocation points. From the field values at 
these points it is also possible to obtain the coefficients 



of a spectral expansion. As in [3J] and [45|] we use stan- 
dard spherical coordinates with Chcbyshcv polynomials 
in the radial direction and Fourier expansions in both 
angles. Within this method it is straightforward to com- 
pute spatial derivatives. To obtain the results described 
below the SGRID code uses at most 3 x 53 + 2 x 161 = 481 
collocation points in the radial and only 64 x 48 in the 
angular directions. This small number of points makes it 
so efficient that it can run on a single PC or laptop. 
As in 13411 we introduce an extra variable 



1 



n = — (d tV j- pd^) 



(46) 



in order to obtain a system of equations that is first order 
in time 



pdiil) - an, (47) 
—^[VrGS'n + aj ij dj^)] + aS cS , (48) 



which results from Eqs. (|46|) and (|40|) . For the time in- 
tegration we use a fourth order accurate Runge-Kutta 
scheme. We implement Eq. (|4"5)l in the code using the 
equivalent, specific form 

d t U = ^dill-ag^didjip + arditp 

- g l ° (diip)dja + aKU + aS cS , (49) 

where K is the trace of the extrinsic curvature of a con- 
stant t hypersurface, and T l is given in terms of the 



In 



Christoffel symbols of the 3-metric as T l = 7 jfe r 
ingoing Kerr-Schild coordinates, K and T l are given by 



K 
V 



1 



H 



3M\ 

(l + #)3/2 r y f — J ' 

1 H /3 4M\ x i 
{1 + H) 2 ~r~ \ 2 + ~r~) T' 



(50) 
(51) 



For the time integration we use a fourth order accurate 
Runge-Kutta scheme. As in [34| we find that it is possible 
to evolve this system in a stable manner if we use a single 
spherical domain, which extends from some inner radius 
i?i n (chosen to be within the black hole horizon) to a max- 
imum radius i? ut- In this case one needs no boundary 
conditions at -R; n since all modes are going into the hole 
there and are thus leaving the numerical domain. At i? ut 
we have both ingoing and outgoing modes. We impose 
conditions only on ingoing modes and demand that they 
vanish. However, since we need more resolution near the 
particle it is advantageous to introduce several adjacent 
spherical domains. In that case one also needs boundary 
conditions to transfer modes between adjacent domains. 
We were not able to find inter-domain boundary condi- 
tions with which we could stably evolve the system (|49[) . 
For this reason we introduce the three additional fields 



4>i = diyb, 



(52) 
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and we evolve the system: 

d t ip = Fdiip-all 

d t U = PdiTL-ag^difa+aPfa 

— g^fadja + aKH + aS e s 
d t fa = J 1 0,o, ■ o,i),.i l <hV,1I \U),<>. (53) 

Note that this system is now first order in both space 
and time and it can be stably evolved using the methods 
detailed below. Also notice that we evolve the Cartesian 
components of all fields. Due to the introduction of the 
additional fields n and fa our evolution system is now 
subject to the constraints 

fa = Btf. (54) 

1. Characteristic modes 
The characteristic modes of the system ([53")) are [47j 

W ± = n±f^ 

w' 1 ' = V- (55) 

For our shell boundaries £j is a spatial outward-pointing 
unit vector. The fields w® and have velocity — j3 % , 
while w have velocity — j3 % ± cv£\ 

2. Domain setup, boundary conditions and initial data 

We typically use 4 adjacent spherical shells as our 
numerical domains. The innermost shell extends from 
R in = 1.9M to r Q = 10M. The next two inter domain 
boundaries are at 18. 1M and 27. 5M. The outermost 
shell extends from 27.5M to R out = 210AZ . The outer- 
most shell always has 161 collocation points in the radial 
direction. The inner shells all have the same number of 
points. We vary their number between 29 and 53. For 
simulations that last longer than about 390M we have 
observed that reflections from the outer boundary can 
reach the particle and introduce errors in the self-force. 
For this reason we have also performed simulations where 
we add an additional outer shell with 161 radial points 
that extends from 210M to R' out = 400M. As we can see 
in Fig. [3j we can now evolve to at least 600M without 
spurious boundary effects. For our simulations we have 
used the Window function in Eq. (TT7|) . 

As mentioned above we do not impose any boundary 
conditions at R ln . At i? ut we impose boundary condi- 
tions in the following way. First we compute d t w + from 
the fields dtip, dtU and dt<fii at the boundary. Then we 
impose the conditions 

d t w~ = — n/r 

dtw^ = P' l (t)i - all. (56) 



3.80x10 




600 



Figure 3: The broken line shows Ft for an outer boundary 
located at R ou t = 210M. We can see reflections arriving at 
the location of the particle at around 390M. If the outer 
boundary is moved out to i?out = 400M (solid line) no such 
effects can be observed during an evolution time of 600Af. 



on the ingoing modes. Finally we recompute dtip, 9*11 
and dt4>i from dtiv^, dtiv® and dtw^ . The motivation 
for the outer boundary conditions in Eq. ([ST?]) is as fol- 
lows. The first equation is equivalent to assuming the 
Sommcrfeld condition ip = f(t — r)/r for some unknown 
function /. The other two conditions are derived from 
the constraints in Eq. ([M]) and can thus be considered 
constraint preserving. 

For the inter domain boundaries we simply compute 
d t w ± , d t w° and d t w^ from d t ip, d t H and d t <j>i at the 
boundary in each domain. On the left side of the bound- 
ary we then set the values of the left going modes d t w~ , 
dtw® and d t w^ equal to the values just computed on 
the right side of the boundary. On the right side of the 
boundary we set dtw + equal to the value computed on 
the left side. This algorithm simply transfers all modes 
in the direction in which they propagate. 

As initial data we simply use ip = II = =0. 



3. Spectral filters 

In order to obtain a stable evolution we apply a filter 
algorithm in the angular directions after each evolution 
step. As in [45| we project our double Fourier expan- 
sion onto Spherical Harmonics. After setting the highest 
I mode in ip and II to zero we recompute all fields at 
the collocation points. This filter algorithm removes all 
unphysical modes and also ensures that ip and II always 
have one less than mode than fa. 



12 



4- Noise reduction 

If we compute the coefficients in a Fourier series expan- 
sion of the effective source for a particle moving along a 
circular orbit we expect them to be of the form 



h m {t) = h m {0)e l 



(57) 



where h m {0) are the coefficients at time t = 0, m is the 
mode number and f2 is the orbital angular velocity. How- 
ever, in the SGRID code we use discrete Fourier trans- 
forms instead of Fourier series, so that the resulting co- 
efficients have a more complicated time dependence for 
any finite resolution. 

The collocation points in the SGRID code are fixed. 
This means that the moving particle periodically ap- 
proaches grid points. Thus for any given resolution, the 
discrete Fourier coefficients of the effective source will 
show a modulation (in addition to the expected phase 
factor) on the timescale it takes to move from one grid 
point to the next. This modulation is a source of ex- 
tra noise. In our simulations we have removed this extra 
noise by the following procedure. We simply compute the 
coefficients h m (0) once and for all at t = 0. For any later 
time we evaluate the source by taking the inverse dis- 
crete Fourier transform of h m (0)e lmnt , so that we avoid 
any extra modulation or noise. 



5. Convergence 

As the source S e g is C° at the particle, we expect that 
ip is C 2 and <pi is C 1 there. This implies that with our 
spectral code ip is expected to be fourth order conver- 
gent at the particle. This expectation is confirmed by 
the results presented in Fig. |U The solid line shows the 
difference in ip between a low and medium angular res- 
olution run, while the broken line shows the difference 
between the medium and high resolution run scaled by a 
factor of s = 3.21 chosen such that the two lines coincide. 
This factor is related to the order of convergence O by 



(1/4. 



(l/N med ) c 



(1/N med )0 - (1/N hi )° ■ 
For an order of convergence of O = 4 we would have 



(58) 



(l/64) 4 -(l/80) 4 
(l/80) 4 -(l/96) 4 



(59) 



The scale factor of s = 3.21 thus corresponds to conver- 
gence of an order between 4 and 5. 



VI. SELF-FORCE AND ENERGY FLUX 

In this section we present and then comment on our re- 
sults. We compute F t and F r , the two non-trivial compo- 
nents of the self-force for a scalar charge in a circular orbit 
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Figure 4: Differences in %p at the particle location for runs with 
different angular resolutions given by Ng. All three cases are 
for iV r = 29 and N<p = 3iVg/4. If we scale the difference be- 
tween medium and high resolutions by 3.21, the two curves 
coincide. This corresponds to convergence of an order be- 
tween 4 and 5. 



around Schwarzschild, and show consistency between the 
results from the two codes. 

Using Eqs. (j2"3"|) and we also compute the scalar 
energy flux across the event horizon and some finite outer 
boundary, referring to this outer boundary as the extrac- 
tion radius. For ease of comparison, these fluxes are ex- 
pressed as the ^-component of the self-force, based on the 
relation given by Eq. (f2"2f . 

A representative summary of the accuracies we 
achieved is presented in Tabic Q] below. 



A. The Dissipative Piece, Ft 

The mode-sum of the t-componcnt of the self-force, 
dt4> K , for the case of a scalar charge in a circular orbit 
in Schwarzschild is known to converge exponentially in 
I, and is thus typically calculated extremely accurately. 
Despite the divergence in i[i rot at the location of the scalar 
charge, dtip rot is smooth there and requires no regular- 
ization. This arises because the retarded and advanced 
fields for a charge in a circular orbit are related by: 



adv 



Writing ifj let 

V> rot = 



(60) 



(61) 



we see that the time derivative of the second term van- 
ishes. The first term is clearly smooth at the location 
of the charge because it is a solution of the homogeneous 
wave equation. For generic orbits, Eq.[5ni w ill not be true, 
and all components of d a yj lct will need to be regularized. 
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Code Result Error 

T t mb (3.728 - 3.748) x 1CT 5 0.05%-0.6% 

F t sgrid (3.7481 - 3.7487) x 1(T 5 0.05% 
T r mb (1.384 - 1.389) x 10 _b 0.4%-0.8% 

F r sgrid (1.384 - 1.386) x 10~ 5 0.4%-0.5% 





Code 


Result 


Error 


E(R = 


150) mb 


3.773 x 10" 5 


0.6% 


E(R = 


150) sgrid 


3.771 x 10~ 5 


0.6% 


E(R = 


300) mb 


3.761 x 10~ 5 


0.2% 


E(R = 


oo) mb 


3.7502 x 10" 5 


0.0005% 



The results of Fig. [5] are from the multi-block code. 
It plots F t evaluated at the location of the charge as a 
function of time t. Helical symmetry would correspond 
to a horizontal line, and we see that the plot gradually 
approaches this while also getting to the correct self-force 
(based on highly-accurate frequency-domain results in 
the literature). The particle makes about two full or- 
bits (T orb = 2it-s/WJM r* 200M) before F t is reached 
to within 1%. The result improves as initial data effects 
further diminish. 



Table I: Summary of (3+1) results. The top half of the table 
reports the computed components of the self- force for a charge 
in a circular orbit r = 10M. These were extracted around 
time, t=600M. The error is determined by a comparison with 
an accurate frequency-domain calculation [l(|, which reports 
F t = 3.750227 x 10" 5 and F r = 1.378448 x 10" 5 . The bottom 
half of the table reports the computed energy fluxes through 
the event horizon and the two-sphere defined by outer ex- 
traction radius R. The R = oo case is an extrapolation to 
infinite outer extraction radius that was performed on results 
of the multi-block code (as explained in Sec. I VI Cj) . For ease of 
comparison with the local self- force, all energy fluxes are ex- 
pressed as Ft according to Eq. (|22p . "mb" and "sgrid" stand 
for multi-block and sgrid codes, respectively. 



By instantaneously switching on our source at t = 0, 
the early part of the evolution will be contaminated by 
initial data effects. After some time these transient ef- 
fects propagate out of the numerical domain and the sys- 
tem settles down to its helically symmetric end state. 
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Figure 5: F t computed at different angular resolutions of the 
multi-block code. The resolutions are described by the num- 
ber of angular gridpoints per patch, so that 40 x 40 corre- 
sponds to 40 gridpoints in both and <j> directions within each 
patch. The noise in this plot is due to the C° nature of the ef- 
fective source. The frequency of this noise corresponds to the 
particle travel time from one gridpoint to the next, and the 
low-frequency modulation corresponds to the particle travel 
time between patch boundaries. 



On top of this evolution towards the correct self-force 
appears to be some sort of modulated noise. This behav- 
ior is the result of two factors. The high-frequency com- 
ponent is due to the fact that the source is only C° (and 
hence derivatives of the scalar field are only C 1 ) at the lo- 
cation of the particle. The finite differencing scheme em- 
ployed here uses stencils near the particle location that 
enclose this non-smoothness, and this is expected to in- 
troduce some noise. The frequency of this noise corre- 
sponds exactly to the particle travel time from one grid 
point to the next. The low-frequency modulation that 
envelopes the noise has a period of about 50M, which is 
exactly the time between crossings of inter-patch bound- 
aries. This is due to the inflated sphere coordinates used 
within the individual blocks. The angular resolution is 
slightly higher near the edges than in the middle of the 
block. 

We observe that the amplitude of this noise decreases 
with increased angular resolution. At 40 x 40 angular 
resolution, we obtain values for the self-force between 
3.7 x 10~ 5 and 3.75 x 10" 5 i.e. within 1.3% of the fre- 
quency domain value (which we will consider in the fol- 
lowing to be exact). At an angular resolution of 60 x 60 
the amplitude of the noise is smaller by a factor of 2.25 
corresponding to second order convergence and an error 
of about 0.6% of the exact value. 

Both results correspond to a radial resolution of Ar = 
A//10. The inner (excision) boundary was placed at 
Ri n = 1.8M and the outer boundary at i? ou t = 400 A/. 
Modifying the radial resolution (to Ar = M/15) does not 
significantly impact the amplitude of the noise. 

In Fig. [6] we compare the results from the two codes. 
There is good agreement between the two, except that 
the SGRID result has noticeably less noise than the 
multi-block result. This is due to the extra noise re- 
duction performed by the SGRID code, as described in 
Sec. IVB4I For the SGRID result shown here, the num- 
ber of collocation points in the 9- and ^-directions were 
Ng = 64 and = 48, respectively. In the r-direction, 
A" 1 = 53 collocations points were used in each of the 
three inner shells. (The number of collocation points in 
the two outer shells, A° ut , was kept the same in all runs 
at A° ut =161.) 
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Figure 6: Comparing Ft results from the multi-block and 
SGRID codes. The multi-block result was achieved with 
60x60 angular resolution and Ar = M/10 radial resolution, 
as in Fig. For the SGRID result shown here, the number 
of collocation points in the angular directions were Ne = 64 
and iV^ = 48. In the r-direction, JV^" = 53, for the three inner 
shells. The two outer shells were always set to use JV° ut =161 
collocation points. 



B. The Conservative Piece, F r 

The conservative piece of the self-force is really the 
crucial quantity to compute. This is the part of the self- 
force that cannot be inferred from observations far away 
(unlike F t , for example, which can be determined from 
the energy flux by using Eq. [22]) . For the case of circular 
orbits, this conservative piece shows up entirely as the 
r-component, F r . In a mode-sum self- force calculation, 
this would be the quantity whose mode sum converges 
as £~ n , where n > 1 is typically a small number depend- 
ing on the number of regularization parameters one has 
access to. In our approach, calculating F r (where r is 
the Schwarzschild radial coordinate) amounts to taking 
derivatives of the regular field ip R , which corresponds to 
taking simple algebraic combinations of the interpolated 
values of the evolved fields at the location of the charge. 

We present the results from the two codes together in 
Fig. [Jj The data from the multi-block code were com- 
puted using runs at 40 x 40 angular resolution and two 
radial resolutions Ar = M/10.M/15. For the SGRID 
code, we have used data from the same run described in 
Fig.H 

We immediately notice that all the results eventually 
settle on a value slightly offset from the correct one. It 
is worth emphasizing though that for the SGRID data 
and the multi-block data calculated at radial resolution 
of Ar = M/15 the final error is just < 1%. Moreover, as 
can be seen from the two multi-block results, this offset 
converges away with increasing radial resolution. 



Figure 7: Comparing F r results from the multi-block and 
SGRID codes. 

C. Energy Flux 

An important consistency check for our runs is the re- 
lation between the scalar radiation flux and Ft, as given 
by Eq. (|22|) . Figures O and [9] display some results from 
the multi-block and SGRID codes, respectively. 

Figure [5] features results from the 40 x 40 angular res- 
olution run of the multi-block code. In this plot, we dis- 
play the energy fluxes through two different outer extrac- 
tion radii, R — 50M and R = 300M, added to the energy 
flux through the event horizon. The outer boundary of 
the computational domain was at R ou t — 600M for both. 
For easy comparison, the energy fluxes are converted to 
a self- force using Eq. ([2"2"| . Also plotted are the results 
from the local calculation of F t (i.e. computed by sim- 
ply taking the time derivative of the regular field at the 
location of the charge) as a function of time also at the 
40 x 40 resolution. These arc all compared with the cor- 
rect frequency domain result represented by the straight 
line. The flux from the larger extraction radius and the 
direct calculation of F t both show agreement to within 
1%. The energy flux is much smoother than the calcu- 
lated local F t , since it is an integral over a spherical sur- 
face of smooth fields far away from the non-smoothness 
at the particle location. 

In Fig. [9[ we see the corresponding results from the 
SGRID code. These come from the same runs described 
in Fig. [5J The energy flux was calculated using an outer 
extraction radius of R = 150M, and again converted 
to the corresponding Ft. This is juxtaposed with the 
local calculation of F t and the frequency-domain result. 
Again, we observe that except for early-time errors due 
to spurious initial data, the energy flux settles to within 
1% of the correct answer. 

One notable observation is that the energy flux im- 
proves with increasing extraction radius. This is shown 
in Fig. 1101 Knowing this, it is tempting to make the ex- 
traction radius as large as possible. However, how far the 
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Figure 8: E computed at 40 x 40 angular resolution of the 
multi-block code. The energy flux includes the contribution 
through the event horizon and an outer boundary defined by 
R. Shown here are results from two outer extraction radii, 
R — 50M and R = 300 M. Energy fluxes are expressed as F t 
according to Eq. (|22|) . Also plotted is the local calculation of 
Ft at the same resolution. 



3.9e-05 



3.85e-05 



3.8e-05 



| 3.75e-05 



3.7e-05 



3.65e-05 



F t (sgrid) 
Edot (sgrid) 
correct F, 



Figure 10: Dependence of E on various extraction radii. Far- 
ther extraction radii is observed to yield better results. 



radiation reaches 50 M first, and only after an interval of 
time arrive at the next extraction radius at R — lOOAf . 

Using the fact that outgoing null geodesies travel at 
coordinate speed (r — 2M)/(r + 2M) in Kerr-Schild co- 
ordinates, one can integrate and find that the time delay 
between the arrival at various radii are given in Tab. [Til 



Interval Time Delay 



50M - 100M 
50 Af - 150M 
50M - 200M 
50M - 250M 
50M - 300M 



47.0351M 
95. 7095 M 
144.572M 
193.687M 
242.963M 



Table II: Time lags. 
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Figure 9: E with the SGRID code from the same run de- 
scribed in Fig. [6] This made use of an outer extraction radius 
of R = 150M. Also plotted is the result of the local Ft calcu- 
lation. 

extraction radius can be situated is limited by the fact 
that the flux taken at farther radii naturally takes longer 
to equilibrate, since the bad initial data waves will have 
to propagate much farther. 

Instead, one can use results from finite radii to extrap- 
olate the energy flux in the limit of an infinite extraction 
radius. This was done with results from the multi-block 
code, and Fig. QT] shows the outcome. The results from 
six finite extraction radii, from R = 50A/ to 30071/ , were 
used to determine the energy flux in the limit of infinite 
extraction radius. First though, one must account for 
the time shifts in the fluxes. Obviously, emitted scalar 



Shifting the data by these appropriate time delays, we 
assume a form for the flux E{R) at finite extraction ra- 
dius R given by: 

N 

E{R) = E(oc) + E + 0(1/R N+1 ), (62) 

n— 1 

Truncating at TV = 5 the constants C±, . . . , Cjv and E(oo) 
can then be determined from the the six sets of flux data 
at the different extraction radii. The resulting E(oo) 
from this procedure is plotted in Fig. 1111 For reference, 
we also include the frequency domain result for F t ex- 
pressed as a flux with Eq. (|22|) . As expected, the agree- 
ment is significantly improved. Extrapolating to infinite 
extraction radius, the flux matches to ~ 0.0005%. We 
take this result as further validation that our effective 
source is a good C° representation for a point charge 
that would otherwise have been represented with a delta 
function. 
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Figure 11: E computed with the multi-block code and ex- 
trapolated to infinite extraction radius. 



VII. SUMMARY AND DISCUSSION 

For our flux values to achieve the reported accuracies 
is noteworthy. This indicates that our effective source 
Sett is a good representation for point particles, in place 
of delta functions that are difficult to handle on a grid. 
Narrow Gaussian functions centered around the location 
of the particle have previously been employed for this 
task. Reference [48| uses a time-domain calculation of 
the gravitational energy flux for a point mass orbiting 
a Kerr black hole which results in errors ~ 10%. But 
by optimizing the number of grid points used to sam- 
ple the narrow Gaussian one can actually get results of 
~ 1% [49( | . Recent work by Sundararajan et al. [5(| has 
done even better than this, with a novel discrete rep- 
resentation of the delta function. Errors of < 1% have 
consistently been achieved with this technique on time- 
domain codes that solve the Teukolsky-Sasaki-Nakamura 
equation. We note that our flux results with a (3+1) sim- 
ulation are already at a comparable accuracy, albeit only 
for the scalar energy flux. It is difficult to speculate on 
how narrow Gaussians and discrete representations will 
perform in a (3+1) context. 

The main highlight, however, has to be the accuracy of 
our self- force results, which have errors < 1%. Unlike any 
other self-force calculations thus far, these values were 
calculated by merely taking a derivative of the regular 
field at the location of the point charge. These results 
are promising as a first attempt at doing a (3+1) self- 
force calculation. 

The judicious placement of collocation points by the 
SGRID code close to the location of the charge appears to 
enable it to represent the effective source better (and to 
achieve slightly more accurate results) as opposed to the 
uniform grid that the multiblock code uses. This seems 
to suggest that devoting more resources to resolving the 
region around the charge (like what would be done with 



adaptive mesh refinement) is the right strategy. 

Both codes show convergence with respect to increases 
in radial and angular resolution. It is clear that the fi- 
nite differentiability of the source reduces the order of 
convergence of the codes relative to what it would be 
if one had a smooth source. For instance, exponential 
convergence ought to be observed in the SGRID code, 
and the non-smoothness of the source is significant for 
the multi-block finite difference code, where only second 
order convergence is achieved while much higher order 
operators are actually used. Modifications of either code 
aimed at treating sources with limited differentiability 
would be likely to improve the order of convergence. Such 
a modification might take the form of an adjustment to 
a stencil or a spectral function to anticipate the location 
of the charge. 

An important issue has been made apparent by these 
initial results. Since the self-force is a very small quan- 
tity, the effects of imperfect boundary conditions become 
a cause of concern. In both the SGRID and multiblock 
codes the outer boundary conditions were implemented 
in a way that ignored the back scattering off the curva- 
ture of the spacetime outside the computational domain. 
Since the self-force contains a tail effect, any such bound- 
ary condition will, when the boundary comes into causal 
contact with the particle location, affect the calculation 
of the self-force. In practice it will seem like the outer 
boundary partially reflects the outgoing waves. In order 
to avoid such effects we have to place the outer bound- 
aries far enough out, that they remain out of causal con- 
tact with the particle (or the sphere where the energy flux 
is measured) for the duration of the run. This of course 
makes the runs more computationally expensive both in 
terms of memory and cpu usage and limits the number of 
orbits that can be simulated. Accurate self-force analy- 
ses require careful treatment of the boundary conditions. 
One way to do this would be to use the non-local radia- 
tion boundary conditions developed by Lau [o"l| and used 
in practice for calculating the metric perturbations of an 
extreme mass ratio binary with a 1+1D discontinuous 
Galerkin code [52j . 

In summary, with this preliminary study, we have 
demonstrated how it is possible to compute self-forces 
with existing (3+1) codes — in fact one of our implemen- 
tations runs on a laptop! Moreover, it has been shown 
that even in the (3+1) context, the effective source is a 
good smeared-out alternative to standard delta-function 
representations of point sources. The flux resulting from 
the effective source matches that due to a point charge 
with very good accuracy. There do remain some ques- 
tions to be explored, like the benefits of optimizing the 
codes, the reduction of the convergence order due to the 
finite differentiability of the effective source, and the lim- 
itations set by the effects of boundary reflections. As this 
is merely a first cut analysis wc shall leave these issues 
for future work. 

A goal of this project is to raise interest within the 
numerical relativity community in self-force analyses of 
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the EMRI problem. Thus the C++ code which evalu- 
ates the effective source for a delta-function scalar charge 
has been placed in the public domain via a website 
http://www.phys.ufl.edu/~det/effsource . Our ini- 
tial expectation is to extend this work by allowing for a 
generic worldline. Our longer term goal is to have code 
for an effective source which represents a point mass or- 
biting a rotating black hole. At each step as the project 
progresses we will continue to put in the public domain all 
of our code necessary for evaluating the effective source. 
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Appendix A: EFFECT OF A C*° SOURCE ON A 
FINITE DIFFERENCE CODE 

In order to better understand the convergence proper- 
ties of a finite difference code for the scalar wave equation 
with a C° source, we turned to the 1+1 D case with unit 
speed in flat space 

d 2 ip 



+ a"2 =s(t, x ). 



(Al) 



Similarly to the 3+1 case we can introduce the additional 
variables p = dtip and <p = d x ip and rewrite the wave 
equation in first-order in time, first-order in space form 



dtp 

d t <t> 



d x <j> - S(t, x) 

dxP 

/>■ 



(A2) 



If V(0, x) = 0, d t ip(0, x) = and if S(t, x) = for t < 
the solution to Eq. (|A1|) at a given coordinate (t, x) can 
be shown to be given in terms of an integral of the source 
over the domain of dependence, i.e. 



ft pX-\-t—t 

i/>(t,x) = -l/2 / S(t',x') dx' dt'. 

JO Jx-t+t' 



(A3) 



The solutions for dtip(x,t) and d x tp(x,t) arc then given 
by 



dtfl>(t,x) = -1/2 / (S(t',x + t-t') 

+ S(t',x-t + t')) dt' 



(A4) 



and 



d x tp(t, x) = -1/2 / (S(t',x + t-t') 

-S(t',x-t + t')) dt' 



(A5) 



These integrals can be evaluated numerically in Mathc- 
matica to high accuracy for any given source thus yielding 
the exact solution to be compared with an approximate 
finite difference solution. 
A function of the form 



f(t,x) = cxp 



x — at 



tanh(x — at) (A6) 



is negative for x < at and positive for x > at. Forming 



s(t,x) = \\f(t,x)\\-f(t,x) 



(A7) 



thus results in a source that is positive for x < at and 
zero for x > at. This source is then C° at x = at and 
C°° everywhere else. 

Solving the system of equations in Eq. (|A2[) with the 
source in Eq. (|A6|) and Eq. (|A7p using fourth order cen- 
tered finite differencing and fourth order Runge-Kutta 
time integrator we can then use the exact solution in 
Eq. (|A4[) to calculate the error in the numerically eval- 
uated p. For the specific choice of source parameters 
a = \/2/2 and c = 1.3 we calculated the solution for 
3 different spatial resolutions Ax = (0.1,0.05,0.025) 
on the spatial interval x G [—6,6]. The timestep was 
At = Ax/A. 

The scaled errors (for second order convergence) in 
p at t = 3 can be seen in Figure [TSl As can be seen 
the errors are as high frequency as can be allowed given 
the spatial resolution, i.e. the error varies dramatically 
from grid point to grid point. Therefore it is impos- 
sible to talk about pointwise convergence since the er- 
ror at a given gridpoint may be positive at one resolu- 
tion but negative at another. However, the amplitude in 
the error can still be considered second order convergent. 
In fact calculating the discrete L2-norm of the error we 
find that ||e(Ax = 0.1)|| 2 /||e(Aa; = 0.05)j| 2 = 4.21 and 
\\e(Ax = 0.05)|| 2 /||e(Ax = 0.025) |j 2 = 4.12 showing that 
we have global second order convergence in the L2-norm. 
The numerical methods used here are formally fourth or- 
der accurate, but because the source is C° we are limited 
to only second order convergence. 

If instead we use the source 



S(t,x) = -f(t,x), 



(A8) 



18 



0.015 



0.01 



0.005 



1 -0.005 



-0.01 



-0.015 



-0.02 




Figure 12: Scaled errors for second order convergence in p at 
t — 3 with a C° source. 




Figure 13: Scaled errors for fourth order convergence in p at 
t — 3 with a C°° source. 



with the same values for a and c as before we obtain 
the scaled errors shown in Figure 1131 Here the errors 
are smooth and low frequency and the scaled errors 
from different resolutions (scaled for fourth order con- 
vergence) agree very well, i.e. we have pointwise con- 
vergence. For the L2-norms of the errors in this case, 
we get ||e(Ax = 0.1)|| 2 /||e(Aa; = 0.05)|| 2 = 15.74 and 
||e(Ax = 0.05)|| 2 /||e(Ax = 0.025) || 2 = 15.92, clearly 
showing the expected reduction in errors by a factor of 
16 with a doubling of the resolution. 



Appendix B: SELF-FORCE AS BOUNDARY 
INTEGRALS IN KERR-SCHILD COORDINATES 



In this Appendix, we derive the relationship between 
the (Kerr-Schild) time-component of the self-force on the 
scalar charge and the energy flux at spatial infinity and 



across the event horizon. We also derive explicit expres- 
sions for these fluxes in Kerr-Schild coordinates. 



1. How Ft and the energy flux are related for 
charges in circular orbits 

For a scalar charge going in a circular orbit around a 
Schwarzschild black hole, there exists a direct relation- 
ship between the time-component of the self-force on the 
scalar charge and the energy flux at spatial infinity and 
across the event horizon. 

In the absence of external fields, the motion of a scalar 
particle is governed by the self-force acting on it, 



ma b = q{g bc + A c )V c f 



F a . 



(Bl) 



The energy per unit mass (i.e. specific energy) of a par- 
ticle along a geodesic with a four-velocity u b is just 
E = — ti,u b , where t b is the time-translation Killing vec- 
tor of the Schwarzschild spacctimc. The rate of change in 
this specific energy per unit proper time is E = u c V ' c E = 
—u c u b V c tb — t b u c V ' c Ub = —t b ab, since V( c ^) = 0. In 
Kerr-Schild coordinates this is just E — —a t . Evaluating 
this on a particle moving in a circular orbit, Eq. (|B1[) 
gives us, 



E\ 



■^(d t i/j R + u t u b V b ^ R )\ p (B2) 
'I 



in 



(B3) 



where p signifies the location of the particle. The second 
term in the first equality vanishes for a circular orbit, 

u a V a i' R \ p = (u*dt V> R + u*d^ R ) \ P (B4) 

= (dt/dT){d t iP R + nd^ R )\ p (B5) 

= (dt/dT)£^ R \ p (B6) 

= 0, (B7) 

where £ a is the Killing vector associated with the heli- 
cal symmetry, so that the third equality vanishes due to 
Eq. Thus 



F t = -mE = qd t ^ R . 



(B8) 



The time component of the self-force in Kerr-Schild co- 
ordinates is then just the amount of energy lost by the 
particle per unit proper time. 

This energy loss must obviously be related to the scalar 
energy flux. We shall now derive these relationships and 
also the explicit expressions for the energy fluxes in Kerr- 
Schild coordinates. 

The scalar field produced by the charge is determined 
by the field equation, 



6^{x-z{t)) 



dr. 



(B9) 
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Multiply both sides by t a V a ip, integrate over V (which 
we take to be the 4-volume bounded by constant Kerr- 
Schild i-surfaces t = i, and t — tf, the event horizon, 
and time-like hypersurface r = R), and t a is the time- 
translation Killing vector of Schwarzschild, and simplify 
the integral over the delta function to obtain 



/ (i a V Q V) V h V 6 V> d 4 x 
Jv 



= -Anq J (to-Va^lpidt/dr)- 1 dt. (BIO) 

We notice now that the integrand on the left can be ex- 
pressed as 



= 4Trt a g ab V c T bc , 



(Bll) 



where T bc is the stress-energy tensor for the scalar field, 
1 ._„ 1 



<J-ibc 

We then have 
t a V c T ca 



(B12) 



-g d x 



Ftidt/dr)- 1 dt, (B13) 



where, specializing to circular orbits, Ft = qdtip\ P . Here, 
we have exploited the fact that dt4> Tet requires no reg- 
ularization and thus equals dtip R . Since t a is a Killing 



vector, V( c t tt ) = 0, and T ca is symmetric in its indices, 
we have t a V c T ca = V c {t a T ca ). Thus, 



V c (t a T c a )V=gd 4 x 



Ftidt/dr)- 1 dt, (B14) 



t a T rr , dZ L 



dV 



Ftidt/dry 1 dt, 



(B15) 



where c?E c is the directional volume element of the 
boundary dV. The integrand of the left hand side is es- 
sentially the conserved current for the scalar field, t a T a c . 

We recall again that for the case of a scalar charge 
in a perpetual circular orbit of angular velocity f2, there 
exists a helical Killing vector £ a given by 



dx a dt 36 



(B16) 



We break up the left hand side of Eq. (|B15[I into the four 
hypersurface integrals, 



t a T r . a dS c = 



dV 



-l c r 2 dX dn + 

r=2M Jr=R 



r c r 2 dtdfl+ I n c Vh d 3 x 



"Vh d 3 x 



t=u 



(B17) 



r 



l a is the null generator of the event horizon, and the rest cancel each other out. This simply means that energy 
of the hatted quantities are the respective outward unit content in each constant-^ hypersurface is the same. But 
normal vectors to the other hypersurfaces making up 3V. we can show this explicitly. 
A is an arbitrary parameter on the null generators k a of 
the event horizon. From the helical symmetry of the 
problem it is easy to see that the last two integrals just 

I 



Consider the time evolution of the total energy in a 
i-hypersurface, 



dt 



3 



2 t a T ac y/r(r-2M) r 2 dr dil = J — {n c t a T ac ) y/r{r - 2M) r 2 dr dil 

3 



-fi / — (n c t a T ac ) y/r{r - 2AI) r 2 dr dft 
Jt ocp 

. , d 



-n 



where we have used the helical symmetry £^F 



(// " C< ° Tac ^ r: ™) r2 sin0 dr dO^j = 0, (B18) 



' ^tIa ) F = 0' f° r an y F. In other words, time evolution is 



dt 1 * b d<p 
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really just axial rotation. 

Thus, for a circular orbit r = r Q , Eq. (|B15|) becomes 



-l c r 2 dX dn 



--2M 



r c r 2 dt dQ 



r=R 



-™ j' F t dt. 



(B19) 



For convenience, we may set the arbitrary parameter 
A on the horizon to be t. If we then differentiate both 
sides with respect to t, we finally get 



dE 
dt 



=2M 



dE 
~dt 



r=R 



3AI 



F t . 



where 



dE 
~dt 



r=2M 



dE 

~dt 



r=R 



t a T ca (~l c )r 2 dfi, 

--2A1 

t a T ca f c r 2 dfi. 



(B20) 

(B21) 
(B22) 



r=R 



These are the general formulas for the energy flux at 
spatial infinity and the event horizon. In the next sec- 
tion, we write them out explicitly in terms of tp and its 
derivatives. 

Finally, Eqn. (|B20|t can be written in the form, 



dE 
~dt 



+ 



r=2M 



dE 

~dt 



r=R 



dE v 
i—r- 

dt 



(B23) 



where E p — —t a u a is the specific energy of a particle 
moving along a geodesic. This is just a statement of the 
conservation of energy: the energy lost by the charge is 
also the energy flowing through r = 2M and r = R. 



2. Scalar energy flux in Kerr-Schild coordinates 



For convenience we write Eqs. (|B21[) and (|B22|) in 
terms of i\) and its derivatives, in Kerr-Schild coordinates. 
These formulas are essentially the same except for their 
unit normals, where one is null and the other spacclike. 

We first note that in Kerr-Schild coordinates, the 
Schwarzschild metric and its inverse are simply 



2M, , 

(]ab = Vab H K a k b , 

r 



-k a k" 



< a = (l,^), fc"=(l,-- 



(B24) 
(B25) 

(B26) 



where r 2 = x 2 + y 2 + z 2 and ?y ab = diag(-l, 1, 1, 1). 

We begin first with the energy flux through the event 
horizon. The event horizon is essentially a surface of con- 
stant retarded time u — t(s)~ r (S)~ 2Mln {r^/2M — 1), 



where the subscript S means that these are Schwarzschild 
coordinates. In Kerr-Schild coordinates these surfaces of 
constant u are 



t = r + AM In {r/2M - 1) + C, 



(B27) 



where C is just a constant. Any particular surface in this 
family can be defined parametrically by the equations, 

t = A, (B28) 

x = r(A)sin6»cos^, (B29) 

y = r(A) sin 8 sin <f), (B30) 

z = r(A)cos(9, (B31) 

where r(A) is defined implicitly by the relation 

A = r + 4Mln(r/2M-l). (B32) 

With this, the null generator of the surface (which is also 
normal to it) is 



/" 



dx a 



1, 



r-2M\ x 



r + 2M 



(B33) 



With the stress-energy tensor given by Eq. (|B12[) and 
using the expressions given in Eqs. (|B24[) - (IB26[) . a small 
amount of algebra yields 



T ab t a r = i> 2 + 



r-2M\ ■ 

r + 2M) V V 



1 (r-2M 

2 { r + 2M 



d c tpd c ij, 



(B34) 



where the overdot means a derivative with respect to t. 
At r = 2M, the energy flux is then simply just 



dE 
~dt 



(B35) 



= -4M 2 * iP 2 dn. 

r=2M Jr=2M 

The normal one-form associated with the hypersurface 
r = R is £ a = d a r = (0, Xi/r). The corresponding nor- 
malized vector is then 



r - 2M 
This leads to the following 



(2M ( 2M\ x l \ , 

, 1 . (B36) 

\ r \ r J r J 



T ab t a f 



a -6 



r-2M 



2M 



2M 



(B37) 
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Thus, the flux through r = R which is just 

"2M 



dE 
~dt 



R' 



r=R 



R 



R - 1M Jr=R 
2M 



R 



1 



R 



ipd r ip 



(B38) 



Taking the limit r — ► oo, this reduces to the more familiar 
flat spacetime case, 

T ab t a f b = ^d^. (B39) 
And so, we have for the energy flux at spatial infinity, 



dE 
~dt 



= lim R' 

R — * oo 



(B40) 



One of the internal checks we perform is to verify that 
Eq. (|B20[) holds by computing the ^-component of the 
self-force and the fluxes given in equations (|B35j) and 
(IB381) . 



To confront the difficulty of the task, we find it advan- 
tageous to introduce the source field tp (r) , where 



for r < r a : p S (r) 



2r- 



(3r*-r 2 ) 



for r > r a : p S (r) = q/r. 



(C2) 



The source field p (r) is completely determined by local 
considerations in the neighborhood of the object, and it 
is chosen carefully to be an elementary solution of 



vV 



— 47T/3. 



(C3) 



Sometimes we call tp the singular field to emphasize the 
q/r behavior outside but near the small source. 

The actual scalar field tp &ct for the problem at hand 
is approximately p near the small object, and the nu- 
merical problem may be reformulated in terms of the 
remainder field 



<^ R = p &ct - p s 



(C4) 



which is then a solution of 



Appendix C: A TOY ILLUSTRATION OF OUR 
METHOD 



Partial differential equations with two dramatically dif- 
ferent length scales are difficult to solve with numerical 
analysis. Consider the example of a scalar field ip of a 
spherical object at rest, centered at r = and with a 
small radius r a . And, the small object has a scalar charge 
density p(r), with p(r) being constant for r < r a and p(r) 
being zero for r > r a . The small object is inside a much 
larger odd-shaped box, and ip = on the surface of the 
box is the Dirichlct boundary condition for tp. For sim- 
plicity assume that spacetime is flat and, with the object 
at rest, there is no radiation and the field equation is 
elliptic. Then 



2,„R 



VV S - 471-p = 0, 



(C5) 



vV 



-47T/9 



(CI) 



where V is the usual three-dimensional flat space gradi- 
ent operator. 

The challenge is to numerically determine the actual 
field tp act as a function of f everywhere inside the box, 
subject to the field equation (|C1|) and the boundary con- 
dition, and then to find the total force on the small object 
which results from its interaction with tp act . 

On the one hand a very fine numerical grid is neces- 
sary to resolve p in and around the object particularly for 
obtaining the force acting on the object. On the other 
hand, a coarse grid would suit the boundary condition 
while speeding up the numerical computation. If the ra- 
tio of length scales is many orders of magnitude, or if 
the small object is represented by a delta function then 
adaptive mesh methods are unlikely to be adequate to 
resolve the small object while using modest resources. 



where the second equality follows from Eq. (|C3[) . The <^ R 
is thus a source free solution of the field equation, and we 
sometimes call it the regular field because the singular 
tp s is removed from the actual field p act in Eq. (|C4[) . 
And if i/> R is determined then simply adding it to the 
analytically known p s provides p &ct . 

A drawback to this formulation might be that the 
boundary condition requires that p K = —p s on the 
boundary of the box. This is likely to be more difficult 
to impose than the original boundary condition. 

A variation of this approach resurrects the original 
boundary condition. We introduce a window function 
W(r) that obeys three properties: 

A. W(r) = 1 in a neighborhood which includes the 
entire source p(r), that is all r < r a . but the neigh- 
borhood might be larger. 

B. W(r) — for r greater than some value r w which 
is not very small. 

C. W(r) has no structure on the small length scale r D 

Then we modify the source field to be tp s = W(r)tp . 
Now the field equation for the regular field (p K becomes 



VV R = -VV S - 4vrp = S(r), 



(C6) 



and this defines a source S(r) that is zero throughout 
the small object, is zero at the boundary and shows no 
variation on the small length scale r a . So, as long as 
W(r) is smooth enough, the boundary condition for <^ R 
is the natural condition that <^ R = 0. 
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In terms of the original source field ip , the source for 
<y5 R is 

S*(Y) = -V 2 (Wip S ) -47T/9 

= -<p s V 2 W - 2W ■ V^ s - V - 4tt/9 
= -<p s V 2 W - 2W ■ Vip s , (C7) 

where the third equality follows from Eq. (|C3|) and from 
property A of the window function. 

In the formulation based upon Eq. (|C6|) . the small 
length scale has been completely removed from the prob- 
lem. The field (p K ought to be relatively easy to evaluate, 
and then the actual field ip act = c/? R + <p s is trivial to de- 
termine. 

This formulation has the bonus that it simplifies the 
calculation of the force on the object from the field. The 
net force is an integral over the volume of the object, 

F = J p(r)V^ act d 3 x. (C8) 

In the original formulation of Eq. (|Clj) . the actual field 
ip act in the integral would be dominated by <p s which 
changes dramatically over the length scale of the object, 
and </? R could be lost easily in the noise of the computa- 
tion. The fact that cp s and p are spherically symmetric 
implies that 

J p{r)Vfi s d 3 x = 0. (C9) 

Then the substitution ip act — ■> tp s + (p n in the integral of 
Eq. (|C8[) leads to the conclusion that 

F = I p(r)V^ R d 3 x. (CIO) 



Thus the force acting on the object depends only upon 
the field ^ R . 

In addition, the field i/? R does not change significantly 
over a small length scale, so if the object is extremely 
small or even a delta function source then an accurate 
approximation to the force is 

F = qV£% =0 . (Cll) 

This redefinition of the problem at hand is broad 
enough to encompass a suggestion by Barack and Gol- 
bourn [3(| to use a window function that is a step func- 
tion of unity in an inner region containing the small ob- 
ject and zero everywhere outside the region. An imple- 
mentation of this idea involves solving for ip act outside 
the region, with the original boundary condition on the 
box, and solving for (^ R within the inner region, with the 
additional boundary conditions that the value and the 
normal derivative of v? R + (p match those of y act on the 
boundary of the inner region. 

If the source p(r) is replaced by a delta function 8{r) 
then if W(r = 0) = 1 and W(r) is C°° with all derivatives 
of W(r) vanishing at r = 0, then ip s = q/r for all r and 
the source S(r) in Eq. (|C7|) is C°°. However, if the nth 
derivative of W is not zero at r = 0, then the source is 
only C"~ 4 . Or, if for some reason the exact expression for 
ip s is not known, and only an expansion is available, then 
again the source may be of only limited differentiability. 

In applications of this approach to problems in curved 
spacetime, the singular field ip s is rarely known exactly. 
This limits the differentiability of the source of Eq. (|C7j) 
which, in turn, limits the differentiability of the remain- 
der </j r at the particle. 
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